function weightToReturn = calculate_metabolites_weight(allMetabolitesInPathway,currentModel, pathwayStatistics,listOfPoolMetabolites)
    all_metabolites_weight = 0;
    
    for i = 1:length(allMetabolitesInPathway(:,1))
        weight = calculate_metabolites_occurence_weight(allMetabolitesInPathway(i,:), currentModel);
        all_metabolites_weight = all_metabolites_weight + weight * 1;
    end
    
    all_metabolites_weight = all_metabolites_weight / length(allMetabolitesInPathway(:,1));
    existing_metabolites_weight = calculate_integral_weight(pathwayStatistics.EXIST_MET_COUNT, 0, pathwayStatistics.MET_COUNT, 1, false);
    pool_metabolites_weight = calculate_integral_weight(length(listOfPoolMetabolites), 0, pathwayStatistics.MET_COUNT, 1, true)
    
    analysisParameters = evalin('base','analysisParameters');
    weightToReturn = ...
        str2num(analysisParameters.allMetabolitesWeight) * round(all_metabolites_weight * 100)/100 +...
        str2num(analysisParameters.existingMetabolitesWeight) * round(existing_metabolites_weight * 100)/100 + ...
        str2num(analysisParameters.poolMetabolitesWeight) * round(pool_metabolites_weight * 100)/100;

    function score = calculate_metabolites_occurence_weight(metaboliteKEGGID, model)
        index = find(ismember(model.metKEGGID, metaboliteKEGGID),1);
        value = nnz(model.S(index,:));
        bound_3 = 2; %rare occurring compounds with
        bound_2 = 10; %medium frequency
        bound_1 = 20; %frequent occurring compounds
        score = 0;
        if value > bound_1
            score = 1;
        elseif bound_2 <= value && value <= bound_1
            score = 1 + (value-bound_1) * (bound_3-1) / (bound_2-bound_1);
        elseif bound_3 <= value && value < bound_2
            score = bound_3 + (value-bound_2) * (bound_2-bound_3) / (bound_3-bound_2);
        elseif 1 <= value && value < bound_3
            score = bound_2 + (value-bound_3) * (bound_1-bound_2) / (1-bound_3);
        end
    end
end